Short title : Chlamydomonas zygote development 1 2 Title : " Gene regulatory networks for the haploid - to - diploid transition of 3 Chlamydomonas reinhardtii " . 4 5

46 47 The sexual cycle of the unicellular Chlamydomonas reinhardtii culminates in the 48 formation of diploid zygotes which differentiate into dormant spores that 49 eventually undergo meiosis. Mating between gametes induces rapid cell wall 50 shedding via the enzyme g-lysin; cell fusion is followed by heterodimerization of 51 sex-specific homeobox transcription factors, GSM1 and GSP1, and initiation of 52 zygote-specific gene expression. To investigate the genetic underpinnings of the 53 zygote developmental pathway, we performed comparative transcriptome 54 analysis of both preand post-fertilization samples. We identified 253 transcripts 55 specifically enriched in early zygotes, 82% of which were not up-regulated in 56 gsp1 null zygotes. We also found that the GSM1/GSP1 heterodimer negatively 57 regulates the vegetative wall program at the post-transcriptional level, enabling 58 prompt transition from vegetative wall to zygotic wall assembly. Annotation of the 59 g-lysin-induced and early zygote genes reveals distinct vegetative and zygotic 60 wall programs, supported by concerted up-regulation of genes encoding cell wall61 modifying enzymes and proteins involved in nucleotide-sugar metabolism. The 62 haploid-to-diploid transition in Chlamydomonas is masterfully controlled by the 63 GSM1/GSP1 heterodimer, translating fertilization and gamete coalescence into a 64 bona fide differentiation program. The fertilization-triggered integration of genes 65 required to make related, but structurally and functionally distinct organelles -66 the vegetative vs zygote cell wall -presents a likely scenario for the evolution of 67 complex developmental gene regulatory networks. 68 69 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. Introduction 70 71 Metazoa, Embryophyta, and Fungi, the so-called called crown groups, 72 independently evolved developmental programs that give rise to complex 73 multicellular organisms. Their explosive morphological diversity is largely driven 74 via modification or reuse of existing gene regulatory networks (GRNs) consisting 75 of transcriptional regulators and signaling pathways (Carroll 2008; Davidson and 76 Erwin 2006; Pires and Dolan 2012; Rudel and Sommer 2003). Therefore, inquiry 77 about the origins of the crown groups has been focused on analyzing the 78 evolutionary precursors of the GRNs that underlie their key developmental 79 strategies. 80 81 Large-scale comparative genomics studies and more focused analyses of 82 ancestral lineages closely related to the Metazoa (e.g. choanoflagellates) and 83 Embryophyta (e.g. charophytes) have indeed reported that many GRN 84 components predate the origins of the crown groups (Fairclough et al. 2013; Hori 85 et al. 2014; King et al. 2008; Mendoza et al. 2013; Wickett et al. 2014; Worden et 86 al. 2009). These findings encourage inquiry into the roles of pan-eukaryotic GRN 87 modules in the differentiation repertoires of the unicellular eukaryotes that share 88 common ancestry with crown-group organisms. Research on the unicellular 89 green alga Chlamydomonas reinhardtii has been particularly useful for gaining 90 insights into the origins of basal GRNs of the Embryophytes, since large numbers 91 of gene families that are not found in fungi or animals are shared within 92 Viridiplantae (Merchant et al. 2007; Worden et al. 2009). 93 94 Given that the most recent eukaryotic common ancestor engaged in sexual 95 reproduction (Goodenough and Heitman 2014), an obvious focus for these 96 inquiries has been the strategies and proteins involved with gametic 97 differentiation, mate recognition, zygote/spore formation, and meiosis. 98 Particularly fruitful have been studies of these processes in unicellular yeasts, 99 where transcriptional and signaling cascades are found to be conserved in the 100 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 2 sexual reproduction of mushrooms and beyond (reviewed in Honigberg and 101 Purnapatre 2003; Neiman 2011; Mata et al. 2002). Sexual development of 102 Chlamydomonas, detailed below, uses several such strategies and proteins. As 103 an example, both Chlamydomonas and yeast gametes contain a transcription 104 factor that forms a heterodimer in the diploid zygote that is involved in diploid 105 development (Goutte and Johnson, 1988; Lee et al. 2008). Furthermore, the 106 Chlamydomonas heterodimeric transcription factors GSM1/GSP1 in the zygote 107 are structural and functional homologs to the KNOX/BELL homeobox 108 heterodimers that are also involved in the diploid development of land plants 109 (Sakakibara et al. 2008; 2013; Horst et al., 2015). Such deeply rooted 110 conservation of sexual GRNs suggests their fundamental importance for complex 111 multicellular evolution. 112 113 The proteins involved in Chlamydomonas sexual development also exemplify the 114 deep ancestry of sexual processes (Speijer et al., 2015). GEX1, for example, has 115 been identified as a gamete-expressed nuclear envelope fusion protein that is 116 required for nuclear fusion in protists, fungi, plants, and many vertebrates (except 117 those, including humans and mice, whose pronuclei do not fuse in the zygote) 118 (Ning et al., 2013). And the primordial gamete fusogen HAP2 likely was present 119 at the origins of sexual reproduction in eukaryotes (Liu et al., 2008). Recently, 120 HAP2 was shown to be a class II fusion protein, in the same family as dengue 121 and Zika virus fusion proteins (Fedry et al., 2017). 122 123 Molecular understanding of sexual development of Chlamydomonas has been 124 driven by molecular genetics studies (reviewed in Goodenough et al. 2007). The 125 Chlamydomonas sexual cycle is invoked by nitrogen starvation, and entails 126 expression of sex-specific proteins involved in preand post-fertilization events 127 during gamete interactions (mating) and zygote formation (Figure 1). Sex-specific 128 (plus and minus) gametogenesis programs are governed by the bi-allelic mating 129 type loci called MTL+ and MTL(Ferris et al. 2002; De Hoff et al. 2013). The 130 FUS1 gene, exclusive to MTL+, encodes a glycoprotein required for the plus 131 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 3 gametes to bind and fuse to minus gametes (Ferris et al. 1996), and the MID 132 gene, exclusive to MTL-, encodes a transcription factor in the RWP-RK family 133 that induces the minus sex program and suppresses the plus sex program (Lin 134 and Goodenough 2007; Ferris and Goodenough 1997). Therefore, both MTL+ 135 and MTLare necessary for a successful progression through the sexual cycle. 136 137 Minus and plus sexual programs employ a series of key players acting in pairs 138 that are expressed in a sex-specific manner: agglutinins (SAD1 and SAG1) on 139 the flagellar membrane (Ferris et al. 2005; Goodenough et al. 1985), fusion140 enabling factors (HAP2 and FUS1) on the plasma membrane (Misamore et al. 141 2003; Liu et al. 2015; 2008), and heterodimeric transcription factors in the 142 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 4 cytosol, GSM1 and GSP1 (Lee et al. 2008; Zhao et al. 2001). SAD1 and SAG1, 143 large hydroxyproline-rich glycoproteins (HRGPs), serve as both ligands and 144 receptors for the initial gamete-recognition event (Ferris et al. 2005), which 145 triggers downstream events collectively called the mating reaction that involves 146 cyclic AMP-mediated signal transduction via a cyclic GMP-activated protein 147 kinase (reviewed in Snell and Goodenough 2009). Responses to cyclic AMP 148 elevation include activation of a periplasmic pro-enzyme called gametolysin (g149 lysin) (Buchanan et al. 1989) that removes cell walls, promotes translocation of 150 additional agglutinins to flagella (Cao et al. 2015), and activates the sex-specific 151 mating structures where HAP2 and FUS1 are localized (Goodenough and 152 Jurivich 1978; Misamore et al. 2003; Liu et al. 2008). During gamete 153 differentiation, the homeobox transcription factors GSM1 and GSP1 are 154 synthesized and stored in the minus and plus cytoplasms; upon gametic cell 155 fusion they heterodimerize and translocate to the nuclei, turning on zygote 156 developmental pathway and preparing fused cells for eventual meiotic 157 germination (Lee et al. 2008). Therefore, the MID-controlled minus and plus 158 programs prepare minus and plus gametes for a series of events, occurring 159 within a few minutes in a step-wise manner, that initiates the diploid phase of the 160 life cycle (Figure 1). 161 162 The transition from gametes to zygotes (the haploid-to-diploid transition) 163 represents a key step in the Chlamydomonas life cycle, analogous to the 164 embryogenesis following fertilization in plants without entailing mitotic growth. 165 Zygote development entails signature events largely completed within the first 12 166 hr, such as zygote-specific wall assembly (Minami and Goodenough 1978; 167 Woessner and Goodenough 1989; Suzuki et al. 2000), selective destruction of 168 chloroplast nucleoids derived from minus gametes and of mitochondrial DNAs 169 from plus gametes (Kuroiwa et al. 1982; Aoyama et al. 2006), flagellar resorption 170 and basal-body disassembly (Cavalier-Smith 1974; Pan and Snell 2005), and 171 nuclear and chloroplastic fusion (Cavalier-Smith 1976). Such a rapid 172 differentiation process would require robust regulatory networks to prevent 173 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 5 promiscuous execution, networks that combine intraand perhaps intercellular 174 cues and coordinate a multitude of events. 175 176 The walls of vegetative and gametic cells serve two core functions. 1) They must 177 protect cells from mechanical and osmotic stresses, and 2) they must be easily 178 removed. Vegetative cells shed the so-called mother cell wall after their mitotic 179 progeny have formed their own cell walls, and upon encountering a gamete of 180 the opposite mating type, gametes must be able to rapidly shed their walls to 181 allow cell fusion. In both vegetative cells and gametes, transcripts for cell wall 182 genes are upregulated when their walls are experimentally released by 183 incubating cells in g-lysin (Kurvari, 1997; Hoffmann and Beck, 2005). One 184 notable difference is that only the vegetative wall needs to support large 185 increases in cell volume, which requires continuous remodeling of the existing 186 wall layers. 187 188 In contrast, the zygote wall has only a single function, which is to provide 189 physical and chemical enclosure that allows Chlamydomonas to survive loss of 190 nutrient sources and extreme environmental conditions, including dessication. 191 The zygote cell wall is unaffected by g-lysin (Schlosser, 1976), but is removed by 192 the zygote hatching enzyme (z-lysin) whose biochemical nature is unknown. 193 194 The zygote wall-specific features include secretion and assembly of zygote195 specific extracellular materials such as (1-3)β-D-glucans into 4-8 complex layers. 196 These are distinct in structure and composition from the vegetative wall (Catt 197 1979; Cavalier-Smith 1976; Grief et al. 1987), which is organized into seven 198 layers entirely made of HRGPs (Roberts et al. 1972; Goodenough and Heuser, 199 1985). In liquid cultures, zygote walls connect neighboring zygotes into extensive 200 sheets called pellicle (Minami and Goodenough 1978; Suzuki et al. 2000). 201 202 Two previous reports document that the haploid-to-diploid transition of C. 203 reinhardtii is initiated by the GSM1/GSP1 heterodimer: 1) ectopic expression of 204 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 6 GSM1/GSP1 heterodimers in vegetative cells can induce early zygote genes, 205 pellicle formation, and ultimately commitment to meiosis (Lee et al. 2008); and 2) 206 zygotes of gsp1 null (biparental31) plus gametes and wild type minus gametes 207 cannot activate early zygotic genes nor degrade minus-derived chloroplast DNA, 208 and they undergo extra fusions (Nishimura et al. 2012). 209 210 The molecular details of the zygote differentiation program have been explored 211 by a transcriptome analysis (Lopez et al., 2015), which reported ~ 600 zygote212 specific genes. However, further details of the regulatory networks for the early213 zygote developmental pathway, and the functional contexts of the zygote-specific 214 genes for zygote differentiation, are largely unexplored. 215 216 The present study analyzes the regulatory network controlling the majority of the 217 early zygote (EZ) transcriptome by considering the mating-induced signaling 218 events and contribution of the GSM1/GSP1 heterodimer. Our results define the 219 predominant role of the GSM1/GSP1 heterodimer, and provide molecular 220 insights concerning how various cellular machineries participate in zygote 221 differentiation. 222 223 224 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 7 Results 225 226 Experimental design for identifying early zygote-specific genes 227 228 To identify which mRNAs derive from bona fide EZ-specific genes, we 229 considered the following in the sample choice for transcriptome sequencing. 230 Firstly, we included haploid and diploid strains in N-replete samples, called 231 TAP1/2 in order to minimize mating-type-dependent or -limited or ploidy232 dependent expression biases. Secondly, genetically diverse strains were used in 233 the two early-zygote samples combining 45 min and 2 hr sets after mating as 234 EZ1/2 in order to minimize strain-dependent biases. We compared our samples 235 with the following published samples representing gene regulatory programs 236 (given in parenthesis) that function prior to zygote development: N-replete and N237 starved samples (program NS) (Miller et al. 2010); plus and minus gametes 238 (program GAM), some being mating-type specific (programs PG and MG) (Ning 239 et al. 2013); and gamete samples following dibutyryl cyclic AMP-treatment 240 (program MR, > 350 genes) or subjected to g-lysin treatment (program gL, > 100 241 genes) (Ning et al. 2013). In the end, the zygote-enriched transcripts that are 242 induced by the NS, GAM, MR or gL programs were excluded to collect bona fide 243 transcripts controlled by the fusion-dependent program (program FD). 244 245 To assess the contribution of the GSM1/GSP1 heterodimer to the FD program, 246 we compared zygote samples in a gsp1 null background (bp31) and its 247 complemented strain (bp31C) (Nishimura et al., 2012). Finally, gamete samples 248 ectopically expressing both GSM1 and GSP1 were included in the above multi249 way comparison as pseudo-zygotes (PZ condition); these cells exhibit zygotic 250 traits such as pellicle formation without undergoing gamete fusion (Lee et al. 251 2008). Hence, the entire dataset included 15 conditions in 11 genetically diverse 252 strains (Table 1 and Supplemental Table S1). 253 254 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 8 Our dataset includes field isolates that are highly polymorphic. Therefore, during 255 the RNA-seq mapping, we allowed 3 mismatches after considering average 256 synonymous polymorphism among the field isolates (nucleotide diversity at 257 synonymous sites, dS=0.0317, Flowers et al. 2015). As a result, mapping rates 258 range between 82.1% and 98.2% of the total reads (Supplemental Table S2). 259 Further details of data processing and analysis are given in Supplemental Note 260 S1. Following read mapping, gene expression values normalized as FPKM 261 (Fragments Per Kilobase transcript per Million reads) were calculated using the 262 V5.3 genome annotation available at Phytozome (http://www.phytozome.net/, 263 Supplemental Table S3). 264 265 Clustering analysis identifies early-zygote specific genes in two distinct 266 pools 267 268 Since defining zygote-specific genes relies on the cross-referencing of all other 269 conditions, we performed un-guided clustering analysis with 13 conditions, 270 excluding two zygote samples of the gsp1 null strain, averaged from 33 271 independent samples following FPKM calculations. Prior to the analysis, 5656 272 (out of 18823) gene models were discarded as Cluster 0 by low expression 273 (having less than 4 FPKM in more than 31 samples), likely representing either 274 pseudo-genes or genes specifically expressed under conditions not included in 275 our dataset. The rest were normalized per model by converting FPKM into 276 percentage values of expression over the sum of average FPKM values under 13 277 conditions, so that condition-specific expression patterns, rather than absolute 278 expression or fold-change, could be highlighted. The resulting percentage values 279 were used to generate 50 clusters of similar expression patterns by Ward's 280 hierarchical clustering algorithm implemented in the Genomics workbench 281 software v4. 282 283 Among the 50 clusters, C33, C43, and C50 show highly up-regulated or nearly 284 exclusive expression in the EZ condition (Box plots, Figure 2A-C), and these are 285 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 9 designated EZ-specific genes. C10 also shows significant up-regulation in EZ, 286 but relatively strong basal expression in all the other conditions as well, and was 287 therefore excluded from the EZ-specific set (Figure 2D). Of special interest is 288 C44, where similar enrichment is found between EZ and g-lysin-treated gametes 289 (PL and ML conditions), in contrast to C24 where g-lysin-induction is strong but 290 EZ-enrichment is absent (Figure 2E-F). Hereafter these two gene clusters are 291 referred to as gL+EZ (C44) and gL-EZ (C24). Detailed FPKM statistics and 292 description of the 50 clusters are available in Supplemental Table S4 and S5. 293 294 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 10 Out of 444 candidate genes from the three EZ-specific clusters (after discarding 295 redundant gene models), we further refined the EZ set by removing genes with 296 less than a 4-fold induction between EZ and the vegetative condition (TAP) (98 297 genes) or with maximum expression of less than 12 FPKM in EZ (54 genes). We 298 further eliminated 39 genes whose zygote-specific expression is only evident in 299 one of the two biological duplicates (EZ1 and EZ2 in Supplemental Table S1). 300 The remaining 253 genes are defined as the core EZ genes (EZ-core, 301 Supplemental Table S6). 302 303 We cross-examined the EZ-core with previously reported early-zygote genes 304 (Supplemental Table S7 for EZ-core annotation). EZY1, ZYS1-4, and ZSP1-2 are 305 found in the EZ-core, while EZY2 is absent due to its low expression. Genes 306 EZY3-23 were identified in a microarray study by Kubo et al. (2008), where EZY 307 genes were divided into two classes based on g-lysin-inducibility. Most of their 308 non-g-lysin class genes (EZY4,7,8,9,15,17,18,19,20,22) are found in our EZ309 core, while their g-lysin-induced genes are mostly found in gL+EZ 310 (EZY11,12,13,14) and in gL-EZ (EZY5,23). EZY3 and 10 in their non-g-lysin 311 class are found in gL+EZ, with clear g-lysin-inducibility both in our data and in 312 Figure 1 of Kubo et al. (2008). EZY6,16, and 21 in their g-lysin class show strong 313 expression elsewhere in our dataset and therefore fall in different clusters 314 (EZY6,16 in C38, and EZY21 in C13). This comparison confirms that our 315 clustering approach collected the majority of EZY genes and precisely 316 distinguished those that are g-lysin-induced. 317 318 Another EZ transcriptome study by Lopez et al. (2015) reported a total of 627 319 zygote-specific gene based on 4-fold up-regulation by RNA-seq, which is more 320 than twice as large as our EZ-core. We found that the difference is mainly due to 321 the inclusion of low-expressed genes (144 models with < 2 FPKM), g-lysin 322 induced genes (60 models), and multiple models for the same gene (30 cases). 323 Of the remaining 393 genes, 82.7% (325 genes) are found among our EZ324 specific clusters and C10 (showing EZ up-regulation with strong basal 325 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 11 expression). Overall, the comparison to the Lopez dataset shows that our 326 clustering method can sort the complex EZ transcriptome into sub-classes and is 327 sensitive enough to collect the majority of differentially expressed genes. 328 329 Expression of EZ-core genes is >94% controlled by the GSP1/GSM1 330 heterodimer 331 332 To learn which genes in the EZ-core are under the control of GSM1/GSP1 333 heterodimers, we analyzed differential gene expression between the zygotes 334 generated from the gsp1 null and its complemented strain (bp31 and bp31C 335 zygotes; Nishimura et al. 2012) using a bioconductor package, DESeq2 (Love et 336 al. 2014). A significant increase in bp31C zygotes compared to bp31 zygotes is 337 interpreted to indicate GSM1/GSP1 dependency. In total, 186 genes show > 4338 fold increased expression (484 genes at > 2-fold increase) and only 11 genes 339 show > 4-fold decreased expression (115 genes at > 2-fold decrease) in bp31C 340 zygotes over bp31 zygotes (Supplemental Table S8, S9 and S10). Those up341 regulated genes are mostly from the EZ-specific clusters (166/186). Of the 253 342 EZ-core genes, 207 (81.8%) show significant > 2-fold up-regulation in bp31C 343 zygotes over bp31 zygotes (FDR <0.05, Supplemental Table S6). In contrast, a > 344 2-fold enrichment of bp31C zygotes was found for only 158 out of 13,185 genes 345 (1.2%) from the non-EZ clusters (excluding C10/33/43/50, Supplemental Table 346 S11). Of the remaining 46 EZ-core genes showing no significant up-regulation in 347 bp31C zygotes, 34 genes are either poorly expressed in bp31C zygotes or found 348 to be g-lysin-inducible, meaning that nearly all (207/219, 94.5%) of the EZ-core 349 genes are dependent on GSP1 for EZ-expression. 350 351 We also examined the differential expression of EZ-core in the pseudozygotes 352 (PZ condition) where the zygote program is activated by ectopic expression of 353 GSM1/GSP1 without gamete fusion. We found that 146 of the 253 EZ-core 354 genes show >2-fold induction (Supplemental Table S6, column Z) over the TAP 355 condition, supporting the dominant role of the GSM1/GSP1 heterodimer in 356 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 12 activating EZ-specific genes. The level of up-regulation in the PZ condition is 357 consistently lower than in the EZ condition, likely due to the fact that the zygotic 358 program is non-synchronously induced in the pseudozygotes following N359 starvation (data not shown). 360 361 To check the reproducibility of our analysis, 10 of the EZ-core genes at high 362 (>200 FPKM), medium (100-200), and low (<100) expression levels 363 (Supplemental Table S12A) were analyzed by quantitative RT-PCR (qRT-PCR) 364 in three biological replicates of 30 min, 1 hr and 2 hr zygote samples and 365 compared to plus and minus gametes (Figure 3A and Supplemental Figure S1). 366 All show early induction at 30 min and continued expression until 2 hr following 367 gametic fusion, in a GSP1-dependent manner. 368 369 The EZ-core is regulated at the transcriptional level via TGAC motif. 370 371 We next asked whether the GSM1/GSP1-dependent gene activation is regulated 372 at the level of transcription given that RNA-seq cannot distinguish between 373 transcriptional and post-transcriptional controls. We constructed reporter 374 transgenes with luciferase coupled to the promoter sequences of two EZ-core 375 genes -g15252(ZSP2A) with very clean EZ-specific expression and 376 Cre06.g256800(SND1A) with additional up-regulation in the activated gametes 377 (MA/PA conditions) -and analyzed their expression in wild type and gsp1 null 378 zygotes (Figure 3B). The ZSP2A promoter drove luciferase activity from 1 hr up 379 to 4 hr after zygote formation (~ 6 fold increase), whereas no activation was 380 found in gsp1 null zygotes. The SND1A promoter showed ~ 100-fold increase in 381 luciferase activity early in wild type zygotes, whereas a delayed but significant 382 ~40-fold increase was found in gsp1 null zygotes. A similar delayed up-regulation 383 of SND1A in gsp1 null zygotes was observed in qRT-PCR analysis (Figure 3A), 384 indicating that the SND1A promoter is controlled by the mating reaction in 385 addition to the GSM1/GSP1 heterodimer. In the EZ-core, we have found 11/253 386 additional mating reaction-induced genes based on the average FPKM ratio of 387 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 13 the MA/PA conditions over the MG/PG conditions (Supplemental Table S6). 388 However, they all show greater expression in the EZ condition, suggesting that 389 their functions are mainly focused on zygote development. 390 391 With the promoter-reporter assay results demonstrating the action of 392 GSM1/GSP1 heterodimer at promoters, we performed in silico analysis to identify 393 putative cis-acting elements in the pool of EZ-specific genes by Amadeus 394 software (Linhart et al. 2008). The analysis identified three 10-base motifs 395 (Supplemental Figure S2A). The most significant hit, "gtGACacGAC," shares the 396 same GACnnGAC consensus with the previously reported ZYgote Responsive 397 cis-acting Element candidates (ZYRE; Lee et al. 2008; Hamaji et al. 2016). The 398 second and third motifs also contain in common a single TGAC sequence, 399 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 14 possibly representing the core motif recognized by the GSM1/GSP1 heterodimer, 400 which is also the motif repeatedly identified for TALE-class homeobox 401 transcription factors with "WFI50N" (also shared by GSM1) (Knoepfler et al. 1997; 402 Krusell et al. 1997). 403 404 Distribution of the three ZYRE motifs was examined to find additional genes 405 under the direct regulation of GSM1/GSP1 heterodimer. Overall, the three motifs 406 are found 6,024 times in 4,468 out of 18,458 genes (24%, Supplemental Table 407 S13). Per cluster distribution shows significant increase in the three EZ-specific 408 clusters C33, C43, and C50 as expected (the highest occurrence being 57% in 409 C43), but also in C10 and C44 (gL+EZ) (Chi-square test, p<0.001, Supplemental 410 Figure S2B, Supplemental Table S14). 411 412 We also asked independently whether C10 (showing strong basal expression 413 throughout the examined samples) and C44 (gL+EZ) genes are directly regulated 414 by the GSM1/GSP1 heterodimer using the differential expression analysis 415 between bp31C and bp31 zygotes. C10 includes 74/475 genes that are 416 significantly up-regulated in bp31C zygotes, and C44 includes 9/162 such up417 regulated genes (> 2-fold, FDR <0.05, Supplemental Table S11). Comparing 418 bp31C/bp31 ratio distribution per our co-expressed gene cluster indicates that 419 both C10 and C44 are indeed significantly higher in the bp31C/bp31 ratio than 420 the average ratio distribution of all the other clusters combined (Wilcoxon-test, 421 p<0.001, Supplemental Table S11). In conclusion, the two statistical analyses 422 suggest that at least some of the C10 and C44 genes are directly up-regulated 423 by the GSM1/GSP1 heterodimer. 424 425 A subset of g-lysin-inducible genes is selectively suppressed by the zygote 426 program 427 428 g-lysin is released during the mating reaction to remove cell walls and allow 429 gametes to fuse. It has long been known that if the supernatant from mating cells 430 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 15 is harvested and applied to vegetative cells or unmated gametes, their walls are 431 also removed, and the cells proceed to assemble new walls via activating genes 432 encoding wall proteins and proteins of unknown functions (Hoffmann and Beck, 433 2005; Kubo et al., 2008; Ning et al., 2013). The walls of vegetative and gametic 434 cells are apparently identical (defined as the VG-wall in contrast to the zygotic 435 wall [Z-wall]; Monk, 1988). Therefore, g-lysin-treated unmated gametes, and 436 newly fused gametes exposed to g-lysin during a natural mating, are a priori 437 expected to up-regulate the same VG-wall transcripts that are up-regulated in g438 lysin-treated vegetative cells. 439 440 As noted earlier, two clusters containing the g-lysin-induced genes present an 441 interesting anomaly: the C44 genes (gL+EZ) exhibit EZ-enrichment, whereas the 442 C24 genes (gL-EZ) do not show any sign of EZ-induction. Moreover, the 443 bp31C/bp31 zygote expression ratio is significantly less than one for the gL-EZ 444 set (Supplemental Figure S3, Supplemental Table S11). This suggests that the 445 gL+EZ subset of g-lysin-induced transcripts persists in early zygotes, whereas 446 the gL-EZ subset does not, either because the transcription of the gL-EZ set is 447 switched off and/or because their transcripts are actively targeted for degradation 448 or decay. 449 450 To explore the differential regulation of gL+EZ and gL-EZ genes, we analyzed 451 several g-lysin-inducible genes in wild-type and bp31 (gsp1 null) zygotes by qRT452 PCR and the promoter-reporter assay used for the EZ-core analysis. Three 453 gL+EZ genes (RHM1, AraGT1, RRA2) and five gL-EZ genes (SEC61G, VSP3, 454 PHC19, GAS28, P4H1) of various expression levels were selected for the 455 analysis (Supplemental Table S12). 456 457 qRT-PCR results show that all three gL+EZ genes, and also SEC61G, were up458 regulated 3~6 fold in gsp1 null zygotes and 5~8 fold in wild-type zygotes (Figure 459 4A). In contrast, the other four gL-EZ genes show significant up-regulation in 460 gsp1 null zygotes but almost no changes or decrease in expression in wild-type 461 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 16 zygotes (Figure 4A). Hence our qRT-PCR results document that the differential 462 regulation of the two gene sets is GSM1/GSP1-dependent. 463 464 The exceptional SEC61G gene actually shows comparable RNA-seq data to our 465 qRT-PCR results, where its expression in bp31C zygotes is 1.5 fold higher than 466 in bp31 zygotes, similar to other gL+EZ genes (Supplemental Table S12). We 467 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 17 manually inspected RNA-seq data of the gL-EZ set (172 genes) and identified 16 468 genes showing comparable expression between the EZ and ML/PL conditions 469 (less than two-fold difference) as potential gL+EZ genes (noted in the column K 470 of Supplemental Table S17), which contains five out of seven ER-residing 471 secretion-related genes including SEC61G in the original gL-EZ set. 472 473 We next sought to learn more about the g-lysin program. We first used our 474 promoter-reporter assay to test if the g-lysin program operates at the 475 transcriptional level. We constructed transgenic strains harboring promoter476 reporter constructs for AraGT1, RHM1, and SEC61G (representing gL+EZ) and 477 for PHC19 (representing gL-EZ); all showed consistent inducible reporter 478 expression when subjected to exogenous g-lysin exposure, indicating that the g479 lysin program regulates both sets via transcriptional activation (Supplemental 480 Figure S4 for the promoter cloning schemes and Supplemental Figure S5 for g481 lysin inducibility test). We then asked how the g-lysin program is affected during 482 zygote development by examining those g-lysin inducible promoters. The RHM1, 483 SEC61G, and PHC19 promoters showed comparable up-regulation during the 484 first four hours after gametic fusion in both wild-type and bp31 (gsp1-null) 485 zygotes (10~15-fold in 3 hr zygotes, Figure 4B). The AraGT1 promoter showed 486 stronger up-regulation in wild-type zygotes, although weak but clear up487 regulation occurred in gsp1 null zygotes (3-fold in 3 hr zygotes, Figure 4B). 488 489 The promoter assay results document that these test genes, and by extension 490 the g-lysin program, operates normally during zygote development and hence is 491 not inhibited by the GSM1/GSP1-dependent program. Therefore, the observed 492 paucity of transcripts produced by gL-EZ genes in wild type zygotes suggests the 493 existence of a post-transcriptional mechanism for their down-regulation, one that 494 does not affect gL+EZ genes. 495 496 Gamete-specific genes are not apparently suppressed by the GSM1/GSP1 497 heterodimer during early zygote development. 498 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 18 499 The zygote developmental program is expected to be accompanied by a 500 negative control of gametic features such as agglutinins and mating structures to 501 prevent additional gametic fusion. We examined expression of known 502 gamete/sex-specific genes in our RNA-seq dataset (Supplemental Table S15). 503 Most known minus-specific genes in C18, including SAD1, GSM1, HAP2, and 504 MTD1 are expressed at comparable levels in gametes and early zygotes. The 505 same trend is true for the plus-specific genes residing in MTL+ such as FUS1 506 and MTA1. This trend indicates that the gamete-specific programs are not 507 apparently suppressed by the zygote program. In contrast, the plus-specific 508 genes residing outside MTL+ such as GSP1 and SAG1 showed rapid down509 regulation as early as 30 min zygotes. Earlier examination of GSP1 and SAG1 510 expressions by qRT-PCR also showed rapid turn-off in early zygotes (Nishimura 511 et al. 2012). Interestingly, SAG1 down-regulation was found equally in both 512 bp31C and bp31 zygotes, suggesting that the down-regulation of the plus513 specific genes may be triggered by cellular fusion rather than the GSM1/GSP1 514 heterodimer. 515 516 Prediction of functional contexts of the early zygote transcriptome 517 518 Our analysis has revealed that EZ-core genes and the gL+EZ set of g-lysin519 induced genes constitute the main body of the early-stage zygote transcriptome 520 dependent on the GSM1/GSP1 heterodimer. Although 203/253 members of the 521 EZ-core were previously identified as zygote-specific genes by Lopez et al (2015) 522 and 30/159 of the gL+EZ as g-lysin induced genes by Ning et al. (2013), these 523 earlier studies offered annotations but not detailed analyses. Therefore, we have 524 gone on to analyze these co-expressed gene clusters, focusing on which 525 functional systems are enriched in the EZ transcriptome and assessing the 526 differences among EZ-core/gL+EZ/gL-EZ gene contents. 527 528 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 19 To analyze molecular details, we performed exhaustive manual annotation as 529 described in Materials and Methods, and categorized the genes into seven major 530 functional groups relevant to the zygote differentiation processes: cell wall 531 synthesis and assembly, organelle-targeted, miscellaneous metabolism including 532 transporters, putative gene regulation, signal transduction, 533 secretion/endomembrane/cytoskeleton system, and ubiquitin system (Detailed 534 annotation results in Supplemental Tables S7, S16 and S17). Figure 5 and 535 Supplemental Tables S18 summarize the functional categories among the three 536 co-regulated gene sets. In the following, we focus on Z-wall assembly where our 537 transcriptome has provided new insights and putative molecular players. 538 539 Molecular components of the zygotic cell wall assembly 540 541 The largest annotated functional category associated with EZ-core/gL+EZ/gL-EZ 542 is cell wall-related genes (46/253, 25/159, and 49/172, respectively), including 1) 543 HRGP-encoding, 2) HRGP-modifying enzymes, and 3) nucleotide-sugar 544 metabolism-related, including sugar-converting enzymes and nucleotide-sugar 545 transporters (Figure 5). 546 547 Cell wall constituents: HRGPs and products of glycosyl transferases 548 Previous studies (Minami and Goodenough 1978; Grief et al. 1987) have 549 documented that the zygotic wall and vegetative wall consist of distinct sets of 550 glycoproteins, primarily HRGPs, recognized by domains dominated by proline (P) 551 and serine (S), often in repeated motifs. While the gL+EZ contains no HRGP552 encoding genes, the EZ-core and gL-EZ include unique sets of HRGP-encoding 553 genes. Of particular interest is the presence of a C-terminal hydrophobic domain 554 in six out of the nine HRGPs in the EZ-core (Supplemental Table S7). Such a 555 hydrophobic tail serves as either a transmembrane domain or a signal for GPI556 like lipid-modification. In silico tools predicting GPI-like lipid modification signals 557 (PredGPI and BIG-PI; Pierleoni et al. 2008; Borner et al. 2003) showed that five 558 of the six candidates likely function as GPI-anchor signal (Supplemental Table 559 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 20 S19). In contrast, none of the 36 HRGPs in gL-EZ displays a hydrophobic C560 terminal tail (Supplemental Table S17). 561 562 A GPI-anchor motif is often associated with plant AGP (arabinogalactan protein)563 type HRGPs but has not yet been experimentally characterized in 564 Chlamydomonas. GPI-type lipid-modification was predicted for a fasciclin-like 565 HRGP, Algal-CAM, in the Volvox vegetative wall (Huber and Sumper 1994), 566 introducing the idea that GPI-anchors or similar structures are likely present in 567 algal cell walls. In support of this idea, most of the GPI-anchor biosynthesis 568 genes are found in the Chlamydomonas genome, of which two enzymes 569 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 21 functioning at rate-limiting steps, PIG-O and PGAP1, are up-regulated in early 570 zygotes (C33) (Supplemental Table S20). 571 572 We found a novel six-Cys-containing homology domain among all the EZ-core 573 HRGPs with C-terminal anchors except Zsp1, which we named MAW 574 (Membrane Associated Wall protein) (Figure 6C for the conserved amino acids; 575 Supplemental Figure S6 for the domain sequence alignment). We searched for 576 additional MAW family members by HMMER3 (Mistry et al. 2013), and identified 577 11 Chlamydomonas members, five members in in Volvox carteri, and seven 578 members in Gonium pectorale (Supplemental Table S19). We could not find 579 MAW homologs in any other published green algal genome sequences. 580 Chlamydomonas MAW proteins have up to three MAW domains, and the 581 majority (10/11) contain Pro-rich domains of PPSPX, PPX, and/or PS repeats, 582 characteristics of HRGP shafts (Figure 6A; Lee et al. 2007). Intriguingly, all MAW 583 members identified are found to possess either one or two trans-membrane 584 domains or a putative GPI-anchor signal at the C-terminus (Supplemental Table 585 S19), indicating that the MAW domain associates with plasma membrane586 anchored cell wall proteins in the Volvocales. 587 588 Phylogenetic analysis of the MAW family reveals six well-supported clades 589 among Volvocalean species (Figure 6B). Clades I and II are distinguished as 590 having a putative GPI-anchor signal at the C-terminus (except CreMAW10). 591 Since all Chlamydomonas genes in Clade I and II are exclusively expressed in 592 early zygotes, we propose that they may serve as key organizers of Z-wall 593 development in Volvocacean algae. Clades IV and V share triple MAW domain 594 configurations, and have diverse expression patterns in Chlamydomonas (Figure 595 6D). 596 The absence of HRGP-encoding genes in the gL+EZ set that contribute both to 597 the VGand Z-walls suggests that VG-wall and Z-wall may be assembled with 598 exclusive sets of materials. We therefore next asked what features are unique to 599 the VG-wall, taking the gL-EZ set as its representative. The gL-EZ includes 36 600 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 22 putative HRGP-encoding genes, 20 of which are pherophorins, the largest known 601 HRGP family in the Volvocales. The Chlamydomonas genome encodes 73 602 pherophorin genes, 28 of which are up-regulated in g-lysin-enriched clusters 603 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 23 (C24/8/34), while 10 are found in gameteand mating reaction-induced clusters 604 (C22/26/41) but none are in EZ-specific clusters (C10/33/43/44/50) 605 (Supplemental Table S21). Hence we propose that pherophorins are hallmark 606 components of the VG-wall. 607 608 HRGP modifying enzymes 609 HRGPs are made from nascent polypeptides that undergo extensive O610 glycosylation on serines and hydroxylated prolines (HYPs) within proline-rich 611 segments during ER-Golgi trafficking. Therefore, production of HRGPs requires 612 prolyl 4-hydroxylases (P4Hs) and various glycosyltransferases. In land plants, 613 distinct glycosyltransferase families are known to build either unbranched O614 glycans found in extensins or complex O-glycans found in AGPs (reviewed in 615 Hijazi et al. 2014). 616 617 Both serineand HYP-glycosyltransferases have been recently identified in 618 plants from studies using plant extensins as substrates: SGT (Peptidyl Serine 619 alpha-Galactosyltransferase) and HPAT (hydroxyproline O620 arabinosyltransferases) belong to the GT8 family and catalyze initial O621 glycosylation on serine or HYP (Saito et al. 2014; Ogawa-Ohnishi et al. 2013), 622 and subsequent arabinosylation is catalyzed by members of two GT77 623 subfamilies, RRA and xyloglucanase-113 (Gille et al. 2009; Egelund et al. 2007). 624 To learn whether HRGP modification capacity or specificity is regulated during 625 wall assembly, we prepared an extensive catalog of the Chlamydomonas 626 homologs for the known HRGP modifying enzymes (Supplemental Table S22). A 627 majority of the HRGP-modification-related genes turned out to be up-regulated 628 by the g-lysinor EZ-program (6/17 P4Hs, and 15/16 HRGP-O-glycosylation 629 enzymes), of which four P4Hs and four O-glycosylation genes are specifically up630 regulated during zygote development. In contrast, we did not find a concerted up631 regulation of N-glycosylation-related genes by the g-lysinnor EZ-program (6/31 632 are up-regulated, Supplemental Table S23). 633 634 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 24 Chlamydomonas HRGPs are also decorated by branched glycans (ex. GP1, 635 Ferris et al. 2001), whose core structure consists of Hyp1-Ara1~3 that are mono636 or di-galatosylated. By contrast, the complex glycans associated with AGPs in 637 land plants consist of a Hyp1-Gal1~3 core, perhaps explaining the absence of any 638 Chlamydomonas homologs to currently identified glycosyltransferases 639 associated with AGP-glycosylation that belong to the GT14, GT29 and GT31 640 families (Liang et al. 2013; Geshi et al. 2013; Nguema-Ona 2014; Ogawa641 Ohnishi and Matsubayashi 2015). 642 643 EZ-specific and g-lysin-induced clusters include many putative 644 glycosyltransferases, mainly of GT47 and GT90 members. Characterized GT47 645 members of Arabidopsis are reported to add sugars to glycans during pectin and 646 hemicellulose biosynthesis (Geshi et al. 2011); however, pectin or hemicellulose 647 has not been found in the Chlamydomonas vegetative wall. Phylogenetic 648 analysis of GT47 family including 47 Chlamydomonas and 39 Arabidopsis 649 members shows a peculiar dichotomy wherein all Chlamydomonas members are 650 related to one Arabidopsis member, At3g57630, which carries an EGF domain, 651 and the remaining Arabidopsis members form separate clades (Supplemental 652 Figure S7). Six GT47 genes are exclusively expressed in early zygotes, and 653 three are found in g-lysin-inducible clusters (C8 and C24). 654 655 GT90 proteins are known as multi-functional xylosyltransferases on mannans 656 when involved in capsule formation in fungi, and as serine-O657 glucosyltransferases when involved in Notch signaling in Metazoa (Acar et al. 658 2008; Reilly et al. 2009), whereas the role of GT90 in plants has not been 659 elucidated. GT90 phylogeny shows great diversity and Chlamydomonas-specific 660 expansion in a number of clades in comparison to Volvox proteins (Supplemental 661 Figure S8). Both EZand g-lysin-induced programs include at least one member 662 of each clade (in red and in blue), suggesting the involvement of GT90 members 663 in both VGand Z-wall assemblies (Supplemental Figure S8). 664 665 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 25 Pellicle formation 666 Pellicle formation in liquid culture begins with wall-to-wall adherence by early667 stage zygotes, followed by extensive build-up of HRGP fibers (Minami and 668 Goodenough 1978; Suzuki et al. 2000). Sulfated polysaccharides or 669 glycoproteins are known to be involved in cell-cell adhesion and to be critical for 670 ECM structural integrity in Metazoans and multicellular Chlorophytes, including 671 Volvox carteri (Misevic et al. 2004; Domozych and Domozych 2014). The EZ672 core contains three genes encoding sulfotransferases with distant homology to 673 Metazoan proteins involved in chondroitin and heparan sulfate biosynthesis; we 674 therefore named them CSR (Carbohydrate Sulfotransferase-Related; 675 Supplemental Table S22E). The CSR phylogeny shows seven clades conserved 676 in Volvocacean algae, and two of the three zygote-specific clades are specifically 677 expanded in the multicellular Volvox, possibly participating in cell-cell adhesion 678 (Supplemental Figure S9). 679 680 The zygote wall components become cross-linked into a matrix that resists most 681 solubilizing agents (e.g. SDS), and evidence for H2O2-mediated crosslinking and 682 isodityrosine (IdT)-linkages in the Chlamydomonas vegetative and zygote wall 683 has been presented in Waffenschmidt et al. (1993), where Tyr-Gly-Gly (YGG) 684 was proposed as a IdT motif. In the EZ-core, we found many YGG motifs in 685 MAW10 with a C-terminal transmembrane domain. A member of the candidate 686 protein family involved in H2O2-mediated cross-linking, glyoxal/galactose oxidase 687 (GOX), is among the most highly expressed genes in the EZ-core. This family 688 has been implicated in cross-linking of plant cell wall constituents such as pectin 689 (K. Sôla, E. Gilchrist; M-C Ralet, C. S, Mansfield and G. Haughn, personal 690 communication) and in direct defense against pathogens via its expected 691 catalysis of H2O2 generation (Zhao et al. 2013). GOX genes are present in 692 Chlamydomonas as a large gene family of 19 members, five of which are up693 regulated by the g-lysinor EZ-program, suggesting their general involvement in 694 cell wall assembly. An interesting feature unique to the Volvocacean GOX is the 695 presence of N-terminal Pro-rich domains (15/19) with repeat patterns such as 696 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 26 PPSPX, PS, or PX, indicating that these GOX proteins are likely part of the cell 697 wall network (Supplemental Table S22F). Phylogenetic analysis of the GOX 698 family shows a dichotomy similar to the GT47 phylogeny, where Arabidopsis and 699 Chlamydomonas have undergone independent expansion (Supplemental Figure 700 S10). 701 702 Support of metabolic demands of glycosylation precursors 703 To support increased HRGP-secretion for zygote wall formation and reported 704 extracellular callose accumulation (Catt 1979; Grief et al. 1987), the cytosolic 705 hexose pool would be in great demand during zygote wall assembly. Hexoses 706 must be activated via UDP or GMP nucleotidylation, and the resultant nucleotide 707 sugars then feed glycosyltransferases (Seifert 2004). A total of 20 genes involved 708 in nucleotide-sugar biosynthesis and 38 putative sugar transporters, mostly in 709 triose-phosphate transporter/nucleotide-sugar transporter (TPT/NST) family, are 710 found in the Chlamydomonas genome (Supplemental Tables S24-S25), allowing 711 us to propose a model for the flow of sugars for glycosylation reactions (Figure 712 7). A notable difference between the Chlamydomonas model and its plant 713 counterpart is the absence of GMD/GER and AXS to produce GDP-L-Fuc and 714 UDP-D-Api, since homologs of these genes have not been identified 715 (Supplemental Table S24). According to our model, the sugar flow via UDP-D716 Glc up to UDP-L-Ara would become highly activated during early zygote 717 development by up-regulation of both NT-sugar transporters (13/34) at non718 plastidic locations and sugar-converting enzymes for UDP-L-Ara (red arrows in 719 Figure 7). Of the NST clades, we noticed size expansion in the clades D, E, and 720 H in both Chlamydomonas and Arabidopsis, where 12 genes (out of 22) are 721 found to be up-regulated in the EZ condition and four genes are up-regulated by 722 g-lysin treatment, suggesting that high demands of glycosylation for zygote wall 723 formation may have driven the expansion of NST family in Chlamydomonas 724 (Supplemental Table S25, Supplemental Figure S11). 725 726 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 27 In addition to the flow of hexose into HRGPs via glycosylation in the Golgi, 727 Chlamydomonas zygotes have another major sink for the hexose pool, callose 728 deposition, contributing up to 30% of all extractable sugars in the zygotic wall 729 (Catt 1979). Two callose synthases and one callose synthase-like gene (GT48) 730 that are exclusively expressed in early zygotes are likely enzymes for such a 731 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 28 huge extracellular build-up (Supplemental Table S26). So far, only trace amounts 732 of glucose have been detected in the Chlamydomonas VG-wall (Catt 1979). 733 734 735 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 29 Discussion 736 737 Our transcriptome analysis comparing wild type and gsp1 null zygotes shows 738 that the EZ transcriptome is mostly, if not entirely, under the regulation of the 739 GSM1/GSP1 homeobox heterodimer, positively at the transcriptional level for the 740 zygote-specific genes, and negatively at the post-transcriptional level for the VG741 wall specific genes. The combination of the two regulatory modules achieves 742 differential regulation of the EZ-core/gL+EZ/gL-EZ gene sets, whose annotation 743 points to molecular differences of the cell-type specific wall assemblies as a 744 hallmark of cellular differentiation. 745 746 Non-guided clustering reveals crosstalk between the EZ and g-lysin-induced 747 programs. 748 Our study design explicitly considers strain variability and multi-way comparison 749 by hierarchical clustering of 13 different conditions and diverse genetic 750 backgrounds. Since zygotes are naturally of mixed genetic make-ups, their study 751 needs to consider genetic variation that may bias the analysis. Therefore, we 752 aimed to identify reliable and robust sets of genes (such as the EZ-core) that 753 warrant more thorough examination. Non-guided clustering provided 754 opportunities to discover expected and unexpected crosstalk among the five 755 transcriptional programs: plus and minus-specific gametogenesis, mating 756 reaction-induced, g-lysin-induced, and EZ-specific. Such a blinded approach has 757 documented that the EZ transcriptome includes a major contribution from the g758 lysin-induced program but not the other mating-associated programs. 759 760 Regulatory network for the assembly of two wall types. 761 The VGand Z-wall assembly programs expectedly share factors involved in 762 HRGP-modification and secretory routes for HRGPs, since both walls consist 763 mainly of HRGPs. Such shared factors are clearly represented by the gL+EZ 764 genes that are expressed equally well during the VGand Z-wall assembly. 765 Following our prediction, the gL+EZ cluster includes 42 ER/secretion-related 766 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 30 genes and 26 HRGP-modifyingand nucleotide-sugar metabolic enzyme767 encoding genes, but no obvious HRGP-encoding genes (Figure 5, Supplemental 768 Table S16). 769 770 Our study defines the EZ-core and gL-EZ genes that are exclusively expressed 771 during the Z-wall or VG-wall assembly respectively, reflecting the fact that the VG 772 wall and the Z wall are very different in their organization and properties. 773 Accordingly, the EZ-core and gL-EZ genes include distinct sets of HRGP774 encoding genes -for example, GPI-anchored HRGPs for the EZ-core and 775 pherophorins for the gL-EZ -indicating that they are indeed the Zand VG-wall776 specific cohorts. 777 778 Interaction between EZand g-lysin programs is highly selective, affecting about 779 half of the g-lysin-induced genes. Follow-up analysis of three g-lysin-inducible 780 promoters provided compelling evidence for persistent expression from these 781 promoters in early zygotes; we also document a subsequent post-transcriptional 782 repression of > 150 genes that is GSM1/GSP1-dependent. Hence the early 783 zygote appears to be programmed to discard the transcripts encoding VG-wall 784 structural proteins that might interfere with the Z-wall assembly, but to maintain 785 production of enzymes for glycosylation and other modifications for cell wall 786 proteins employed in both VG-wall and Z-wall construction. 787 788 The post-transcriptional mechanisms employed by the early zygote to recognize 789 and discard VG-wall-specific transcripts remain to be investigated. One 790 possibility is suggested by the presence of Tudor Staphylococcal Nuclease 791 (TSN) in the gL+EZ set (Supplemental Table S19). The TSN homolog in 792 Arabidopsis has been reported to be involved in regulating the stability of the 793 mRNAs that are targeted to the ER for co-translational translocation (Sundström 794 et al. 2009); more recently, its regulatory action has been found to depend on its 795 nuclease domain and on selective decapping of mRNAs at the stress granule 796 (Gutierrez-Beltran et al. 2015), where TSN operates together with a decapping 797 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved. 31 enzyme like Dcp1 and a 5'-3' exonuclease like Xrn1, proteins that are also 798 conserved in Chlamydomonas. 799 800 Departure from unicellular contexts 801 In flowering plants, the haploid-to-diploid transition is followed by sporophyte 802 development, giving rise to multiple tissue and organ types, whereas the haploid 803 phase is dramatically shortened, entailing only a few mitotic divisions. In the 804 moss Physcomitrella patens, key players involved in the haploid-to-diploid 805 transition include KNOX-family and their putative binding partners BELL-family 806 transcription factors, necessary for spore maturation and suppression of the 807 haploid program, and Clf /Fie as a component of the PRC2 that suppresses the 808 diploid program (Sakakibara et al. 2008; 2013; Mosquna et al., 2009; Okano et 809 al. 2009; Horst et al., 2015). Since GSM1 is the only KNOX-family transcription 810 factor in Chlamydomonas (Lee et al. 2008), it is of great interest to understand 811 how the GSM1/GSP1-dependent program evolved into the KNOX-dependent 812 regulation of land-plant haploid-to-diploid transition. A recent study has shown 813 that the two divergent KNOX subfamilies in land plants function in antagonistic 814 roles: type II promotes differentiation programs whereas type I counteracts type II 815 by promoting the undifferentiated state of stem cells (Furumizu et al. 2015). A 816 plausible scenario, therefore, is that the type II KNOX retains the ancestral 817 differentiation-triggering function of GSM1 and that type I evolved to delay the 818 differentiation program leading to meiosis, thereby allowing the emergence of a 819 multicellular diploid phase. Recently, Frangedakis et al. (2016) reported that the 820 type I and II KNOX already diversified in charophyte algae before the evolution of 821 multicellular diploid generation. The dual-action mechanisms of the GSM1/GSP1 822 heterodimer described in this study, controlling both positive gene expression 823 and post-transcriptional mRNA degradation, may have opened up an 824 evolutionary path for the KNOX-based GRN to become bifurcated by sub825 functionalization following gene duplication. 826 827 Materials and Methods 828 www.plantphysiol.org on July 14, 2017 Published by Downloaded from Copyright © 2017 American Society of Plant Biologists. All rights reserved.